Sampling along reaction coordinates with the Wang-Landau method 
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The multiple range random walk algorithm recently proposed by Wang and Landau [2001, Phys. 
Rev. Lett., 86, 2050-2053] is adapted to the computation of free energy profiles for molecular systems 
along reaction coordinates. More generally, we show how to extract partial averages in various 
statistical ensembles without invoking simulations with constraints, biasing potentials or unknown 
parameters. The method is illustrated on a model 10-dimensional potential energy surface, for which 
analytical results are obtained. It is then applied to the potential of mean force associated with 
the dihedral angle of the butane molecule in gas phase and in carbon tetrachloride solvent. Finally, 
isomerisation in a small rocksalt cluster, (NaF)4, is investigated in the microcanonical ensemble, 
and the results are compared to those of parallel tempering Monte Carlo. 
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I. INTRODUCTION 

The calculation of free energies, or free energy differ- 
ences, is a difficult task for simulation. Free energies 
cannot be conveniently expressed as averages, and be- 
cause of this specific answers have been proposed, in- 
cluding the well known thermodynamic perturbation and 
integration methods (|, 0] , but also biasing potentials || 
and umbrella sampling J3J, [5| , constraint dynamics in the 
so-called "blue-moon" ensemble lq, m, adiabatic switch- 



and the recent adiabatic 
, |ll[ . Given two states A 



ing H and reversible-scaling 
molecular dynamics method 
and B in configuration space, these techniques allows one 
to get a thermodynamic picture of the reaction pathway 
between A and B, but also to estimate transition rates. 
When the reaction coordinate is not well defined, it is 
still possible to calculate such rates using the transition 
path sampling approach |l2|, 

Most of the above methods suffer from important limi- 
tations when it comes to a practical situation. One must 
either find by trial and error a suitable guiding poten- 
tial (the multicanonical ensemble sampling jl4| can help 
in this search), or deal with the specific transformation 
rules required by constraint dynamics, which can often 
be cumbersome in the case of angle coordinates. The adi- 
abatic free energy dynamics method |ll| also makes 
use of variable transformation and particular integrators. 

Recently, Wang and Landau have proposed a very sim- 
ple algorithm to compute the density of states of spin 
systems [jL5j . Their multiple range random walk method 
has since been extended to spin glasses M|, continuous 
fluids (T^j , proteins |Tj| and polymer fluids iQ . A com- 
bination with the kinetic Monte Carlo algorithm has also 
been made by Schulz, Binder, and Muller [^0|. Although 
this method was originally conceived to calculate a mi- 
crocanonical quantity, constant temperature properties 
can be recovered using appropriate Laplace transforma- 
tions. The main interest of the algorithm is that it does 
not require any information a priori, other than a suit- 
able choice of reaction coordinate. Instead it converges 
to the expected solution with a greater accuracy at each 
new iteration. The goal of this paper is to propose a sim- 



ple way to use the Wang-Landau method in the context 
of free energy profiles. As will be shown below, the lim- 
itations mentioned earlier essentially vanish, making the 
method really useful in various statistical ensembles. 

The article is organised as follows. In the next section, 
we give the main ideas of the Wang-Landau algorithm in 
the framework of free energy calculations. In section III, 
we illustrate this method on three examples. Firstly, a 
simple analytical 10-dimensional potential energy surface 
is investigated and the simulation results are compared 
with analytical data. Secondly, the free energy profile 
of the butane molecule along the dihedral angle coordi- 
nate is calculated in gas phase and CCI4 solvent, and 
compared with previous works on the same system pl| . 
Thirdly, we study the cube-ring isomerisation in a small 
rocksalt cluster, namely (NaF)4, in the microcanonical 
ensemble. We briefly summarise and conclude in sec- 
tion 



IV 



II. 



THE WANG-LANDAU METHOD FOR FREE 
ENERGY PROFILES 



We consider a classical N-dimensional system in a given 
statistical ensemble characterized by a equilibrium distri- 
bution p. This distribution is a function of the cartesian 
coordinates R = {a;,-, yi, z{\, and also of external param- 
eters such as temperature or total energy. A reaction 
coordinate A(R) is constructed, and we look for the prob- 
ability density p(Ao) for the reaction coordinate to take 
the particular value Aq: 



KAo) = (<5[A(R) - A ]>. 



(1) 



In this equation, (A) denotes the average value of A in 
the statistical ensemble characterised by p: 



(A) = J A(K)p(R)dR I J p(R)dR. 



(2) 



The original multiple range random walk method of 
Wang and Landau consists of performing a Monte 
Carlo simulation using a Metropolis acceptance proba- 
bility with a dynamical weight g(E), where E(R) is the 
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energy of configuration R: 

acc(R old -» R ncw ) = min 



9(E, 



old) 



(3) 



For each visited state with energy E, g is scaled by a 
constant factor /: g(E) — ► / x g(-E). The scaling factor 
/ is initially quite large (~ 2-2.5), and is periodically re- 
duced toward one. g is initially set to one for all values 
of E. As shown by Wang and Landau G3|, this algo- 
rithm smoothly converges to a flat probability distribu- 
tion when / —* 1, meaning that the function g converges 
to the density of states f2 of the system. 

Since Q(E) and p(A) play similar roles, the same ideas 
can be applied to the problem of free energy profiles along 
reaction coordinates. Let 17(A) be a function initially set 
to one in the whole range of accessible values of A. For 
practical and numerical purposes, it turns out that it is 
more convenient to work with s(A) = 1115(A) instead of g. 
We perform a Monte Carlo simulation using the following 
Metropolis acceptance rule: 



acc(R, 



old 



Rn 



p(R„ 



p(Roid) 

P(Riicw) 



/) ^ ff(Aold) 

X 



.9(A now 

e -*(Anew) 



p(Roid) e 



-s(X 



(•4) 



with A ncw = A(R now ) and A i d = A(R id)- After the 
new configuration is visited, the corresponding s is up- 
dated: s(A) — > s(X) + a. The parameter a = In / is pro- 
gressively decreased to zero at each new iteration. The 
histogram h(X) of visited values of A is considered flat 
when h(X) > e(h(X)} for all A. The parameter e is usu- 
ally taken between 0.5 and 0.99. After a large number of 
MC steps and some iterations, g = e s gets close to the 
probability density p defined in equation (Q). As noted 
by Wand and Landau Jig ], the continuous change in the 
statistical weight in Eqn. (||J) strictly violates detailed 
balance. Actually, detailed balance is satisfied only after 
the function s(A) becomes nearly constant, and within 
an accuracy of order of magnitude a |L5| |. In fact, as is 
well known in statistical mechanics, detailed balance in 
only a sufficient condition for the Monte Carlo chain to 
be Markovian 0, §. 

In the canonical ensemble at temperature T = l/feg/3, 
the equilibrium distribution is given by the Boltzmann 
factor p(R) = exp[— (3V(R)], where V is the potential 
energy of configuration R. Introducing T(A) = — s(X)//3, 
equation (0) can be rewritten as 



acc(R, 



old 



Rn 



min[l,exp(-/3AF)], 



(5) 



with AF — F(H ncw ) — _F(R id), and the Landau free en- 
ergy F(R) = V{R) - r[A(R)]. The function T converges 
to the potential of mean force (PMF) W: 



W(X) = -~lnp(X). 



(6) 



In the microcanonical ensemble at total energy E, the 
equilibrium density is p(R) = [E — ^(R)]^ 2-1 , and s(A) 



plays the role of a Landau entropy. Extra variables can 
be included, as in the isothermal-isobaric NPT ensem- 
ble, or the grand-canonical /iVT ensemble. Expressions 
similar to equation (|^) are then found for the acceptance 
probability, based on their usual forms without reaction 
coordinates [jl], ||. 

There is no other input to this algorithm than the func- 
tion g, initially unknown and set to one. In order to illus- 
trate its efficiency, we have chosen to apply the method 
to a variety of molecular systems. 



III. SIMULATION EXPERIMENTS 
A. Ten-dimensional model potential 

Following Rosso and cowor kers (H],|n), we start to test 
our method on a simple potential energy surface (PES) 
for which exact results can be obtained. The expression 
for V(X.) — V(xi, . . . ,xio) is that of a double-well po- 
tential 



V(X) 



D{x\ 



If 



1 



10 



5 I* 

i=2 



10 

X\ ^ y CLiXi . 
i=2 



(7) 



The parameters are D = 5 and a.; = 1 for all i. The reac- 
tion coordinate is A = x\, and the PMF at temperature 
T is found to be 



P7(A)^(A 2 -l) 2 -^£ 



(8) 



therefore W does not depend on T. We have calculated 
the free energy profile W along the reaction coordinate 
A in the range — 2 < A < 2 using the Wang-Landau 
method, with 20 iterations of 1.2 x 10 6 Monte Carlo steps, 
of which the first 2 x 10 5 were removed for equilibration. 
The histogram in A is made of 1000 bins. To illustrate 
the power of the algorithm, the temperature is chosen as 
only T = 10~ 3 . In figure [j] we have plotted the calculated 
PMF against the analytical result of equation (|^) . The 
two curves are nearly indiscernable, and inspecting the 
error function 8{X) = Wo X act(A) — Wsimui(A) in the inset of 
this figure reveals how accurate the computed free energy 
profile is. In the microcanonical ensemble, when the total 
energy E is close to the saddle energy V(0) = D, the true 
dynamics is strongly slowed down and the crossing rate 
between the two wells drops to zero. The Landau entropy 
s associated with the microcanonical probability density 
can be easily calculated as 



s(A) = (N — 3/2) In 



A 2 



E-B(A 2 -l) 2 -yE fl 



i>i J 



. (9) 



up to an additive constant. The simulation was per- 
formed with the microcanonical distribution p(X) = 
[E — y(X)] 4 , with the same number of MC steps and 
iterations as above, at the total energies E = 20, 10, and 
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FIG. 1: Potential of mean force W(X) of the 10-dimensional 
potential of equation ([j]) at temperature T = 10 -3 , for the 
reaction coordinate A — x\. Inset: absolute error between 
simulation and exact results, AW = Wcxact — Wsimui- 



C/5 




FIG. 2: Landau entropy s(A) of the 10-dimensional poten- 
tial of equation (Q) in the microcanonical ensemble at total 
energies E = 20, E = 10, and E = 5.5, respectively. The sim- 
ulation results are labelled (S), the exact results are labelled 
(E). 



B. Gas phase and solvated butane 

We turn now to a more realistic system and a more 
complex reaction coordinate. The cis-trans isomerisation 
in n-butane has previously been investigated by several 
authors pj ^2| . In particular, Depaepe and cowork- 
ers have used this molecule as a benchmark for free en- 
ergy calculations ^l). The accurate results obtained by 
these authors using umbrella sampling and constrained 
dynamics in the blue-moon ensemble allow a stringent 
test of the present method. 

The butane molecule is modelled with a united-atom 
potential interacting with the solvent molecules via a sim- 
ple Lennard- Jones (LJ) potential. The total potential 
energy V of the system having one solute molecule and 
N solvent molecules is written in cartesian coordinates 
R = {ri, . . . , rAr +4 }, where {n, . . . , r 4 }, label the CH 3 
and CH2 groups of the butane molecule, and r, , i > 4 la- 
bel the solvent molecules. The CH3 and CH2 groups are 
connected by harmonic C-C bonds and C-C-C bending 
forces. A torsion potential between the extremal methyl 
groups consists of a dihedral angle term and a LJ inter- 
action U(,(ri4). The solvent is chosen as carbon tetrachlo- 
ride, which is nonpolar and simple enough to allow for 
a simple united-atom description with the LJ form. The 
expression for V is 

V(R) = Vbutanc(ri, . . . ,r 4 ) + V r so lvent(r5, ■ ■ ■ , ""jV-nO+Vint (R) , 

(10) 

with 

3 k 2 k 

y butanc (n, . . . ,r 4 ) = J2 T (|d 4 |-d*) 2 +E ^(e t -e*) 2 +u 3 cos(3a) 



i=i 



Kolveixt (r 5 , . . . , rAr +4 ) = u s (r 



4<Kj 



Vi n t(R) = ^~]y^ Ubs(n, 
i=l j>4 



(11) 

(12) 
(13) 



In equation (11), |dj| and Qi are the C-C bond lengths 
and C-C-C bending angles, respectively: dj = r i+ \ — r,, 
cos#i = — dj.dj_|_i/|dj|.|d.j+i|. The dihedral angle a is 
defined as zero in the cis conformation: 



cos a 



(di Ad 2 ).(d 2 Ad 3 ) 
Idr Adal.lda Ad 3 | ' 



(14) 



5.5, respectively. The results are represented in figure ||. 
The agreement between the numerical experiment and 
the exact value is again very good, even for the most dif- 
ficult case E — 5.5. Therefore the Wang-Landau method 
seems very efficient for computing free energies and re- 
lated quantities along reaction coordinates, at least in 
such low dimensional model problems. 



All parameters are taken from Ref. [|2l) for sake of 
comparison. We recall them for clarity: d\ — d% = 
1.54 A, d* 2 = 1.52 A, k s = 1882.8 kJ moh 1 A~ 2 , 
9\ = 0% = 1.937 rad, k h = 376.56 kJ mor 1 rad~ 2 , 
1*3 = 6.6944 kJ mol -1 . The LJ parameters for the 
Ubifu) term in equation ( p"T| ) arc e = 0.4184 kJ mol -1 
and a = 3.385 A. For CH 3 , CH 2 , and CC1 4 , they are 
respectively e = 0.8661, 0.4778, and 3.087 kJ mol^ 1 , 
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and a = 3.775, 3.982, and 5.27 A. The parameters 
for mixed interactions are obtained from the Lorentz- 
Berthelot combination rules. 

Monte Carlo simulations were carried out at room tem- 
perature T — 300 K, for the n-butane molecule in gas 
phase and in CCI4 solvent at density 1.614 g cm -3 . As in 
Ref. pH l, 123 carbon tetrachloride molecules were placed 
in a cubic box of side L — 27.124 A, periodic bound- 
ary conditions were implemented through the minimum 
image convention. All LJ interactions were truncated at 
11.75 A. 

The reaction coordinate is the dihedral angle a, and 
all simulations consisted of 20 iterations of 1.2 x 10 6 MC 
cycles each, including 2 x 10 5 initial thermalisation cy- 
cles. 1000 bins were used for the histogram in a. In the 
case of the solvated molecule we have supplemented the 
present algorithm with the preferential sampling scheme 



of Owicki and Scheraga 



as implemented by Mehro- 



tra and coworkers j24j . Each carbon group of the butane 
molecule was chosen for displacement with probability 
uii = 1/8. The probability of choosing a solvent molecule 
i decays with its distance from the centre of mass of the 
solute molecule as Wi = 7/rj^, where 7 is determined 
by the normalisation of the total probabilities. Because 
of this bias, the acceptance probability (pj) is modified 
and becomes [E3f 



acc(R, 



old 



R ncw ) = min 



"/new / /j a 7— 1\ 

1, exp(-/3AF) 

Wold 



(15) 



Preferential sampling is particularly useful in this prob- 
lem, because we focus on a property of the solute 
molecule influenced by the solvent. The probability den- 
sities of the dihedral angle g(a) in gas phase and CCI4 
solvent are plotted in figure |3J. The shape is very sim- 
ilar to the curves reported by Depaepe et al. Q, but 
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FIG. 3: Normalised probability density of the dihedral angle 
in the n-butane molecule at T — 300 K, in carbon tetrachlo- 
ride solvent (solid line) and in gas phase (dashed line). 

quantitative comparison must be made on specific indi- 



cators. We have calculated the fraction of trans isomer, 
Xt = J*i20° 9( a )d a , an d the absolute probability den- 
sities j(a') at the free energy barrier — 120°. The 
results are given in table [i]. The agreement between the 

XT g(^) x lOVdeg" 1 



0.664 ±0.002 
0.659 ± 0.004 



0.670 ±0.011 
0.671 ±0.016 
0.668 ±0.020 



9.42 ± 0.07 
9.49 ± 0.10 



9.69 ±0.19 
9.86 ±0.36 
9.68 ±0.35 



Gas phase 
Present work 
Umbrella sampling" 

CCL4 solvent 
Present work 
Umbrella sampling" 
Blue-moon ensemble" 

a Reference |n| 



TABLE I: Distribution of trans conformation of n-butane at 
equilibrium. \t = J 120 o 9( a )da is the fraction of trans iso- 
mer, g(c^) is the absolute probability density of finding the 
system at the free energy barrier ot 1 = 120°. The errors cor- 
respond to one standard deviation. 

work by Depaepe and coworkers and the present results 
using the Wang-Landau method is quite remarkable, es- 
pecially considering that the numerical cost of the Monte 
Carlo simulation is lower than the molecular dynamics 
simulation of Ref. |2l[] . This quantitative agreement sug- 
gests that the Wang-Landau method is accurate even for 
large molecular systems. 



C. Cube-ring isomerisation in a rocksalt cluster 

The previous results have shown that free energy pro- 
files can be efficiently calculated with the multiple range 
random walk technique. Data of comparable accuracy 
could be obtained with the adiabatic free energy dynam- 
ics of Rosso and coworkers |l(], |llj , which has been espe- 
cially designed for use in the canonical ensemble. We now 
illustrate our method on another problem of molecular 
character, but in the microcanonical ensemble. The small 
(NaF)4 ionic cluster has two stable isomers, a cube (Td 
symmetry, energy E = —1.2317 Hartree), and an octag- 
onal ring (C4U, E = —1.2192 Hartree). The cube=^ring 
isomerisation dynamics has been previously investigated 
in a similar (NaCl)4 cluster by Heidenreich, Oref, and 
Jortner f25"f . This system is made of 8 Na + and F~ ions 
interacting via Coulomb forces and Born-Mayer short- 
range repulsion: 



V(R) = J2 A H exp(-/»yry) ± ^ 



i<3 



(16) 



The parameters are taken from ab initio calculations p6[ : 
A ++ = = A+_ = 41.777203 Hartree, p ++ = = 
= 0.517745 bohr -1 . A suitable reaction coordinate 
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is required to distinguish the two isomers at a given total 
energy. After several attempts, we have found that the 
principal momenta of inertia could achieve this goal. We 
choose the following reaction coordinate: 



A(R) = 5> 4 -r G | 



(17) * 



where rc is the position of the cluster centre of mass. A 
has roughly the value 100 bohr 2 in the cube isomer, and 
about 200 bohr 2 in the ring isomer. The MC simulations 
were still performed for 20 iterations of 1.2 x 10 6 cycles 
each, with 1000 bins in the range 80 < A < 250 for the 
histograms. The probability densities of A in the micro- 
canonical ensembles at total energies E = —1.2, —1.195, 
— 1.19,-1.18, —1.16, and —1.14 Hartree are plotted in 
figure [|. The ring isomer is much more "flexible" than the 
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FIG. 4: Normalised probability density of the reaction coordi- 
nate [equation ([n])] in (NaF)4, in the microcanonical ensem- 
ble at total energies in the range —1.2 < E < —1.4 atomic 
units. The isomers corresponding to each of the two peaks 
are also represented. 

cubic isomer, which is reflected in the relative widths of 
the two peaks in each of these distributions. The fraction 
of cubic isomers continuously decreases in this energy 
range, and this can be quantified using \ — Jgo ° 9WdX. 
In figure [s] we represent the variations of \ versus to- 
tal energy. For comparison, we have also calculated x 
from parallel tempering Monte Carlo |?7| in the micro- 
canonical ensemble pq| . For this, we also recorded the 
probability distributions of A, and then we extracted x- 
The two methods give very similar results, within the er- 
ror bars, and show that the cube and ring isomers are 
in equal ratios at E ~ —1.196 Hartree. From this study, 
it could be also possible, in principle, to estimate tran- 
sition rates that could be compared with actual molecu- 
lar dynamics simulations. However, one must be careful 
that molecular dynamics simulations do not strictly sam- 
ple the microcanonical ensemble. Actually a geometrical 
factor, which depends on the inertia momenta, hence on 




-1.22 



-1.20 



-1.18 -1.16 -1.14 

Total energy (hartree) 

FIG. 5: Fraction of the cube isomer x = J^q° ff(A)dA of 
(NaF)4 in microcanonical equilibrium. The results of the 
Wang-Landau method (full circles) are compared with the 
results of parallel tempering Monte Carlo (empty triangles). 
The dashed lines indicate the equilibrium point where both 
fractions are 50%. 



A, must be incorporated in the equilibrium statistical dis- 
tribution p(R) p9| . In order to compare with MD, the 
present MC simulations should be carried out again with 
this extra factor in the acceptance probabilities. 



IV. CONCLUSION 

The Wang-Landau method is successfully applied to 
the problem of calculating free energy profiles along re- 
action coordinates. This algorithm has several advan- 
tages over most other schemes currently in use. First, 
it does not require any use of constraints, which often 
imply complex variable transformations and the calcula- 
tion of the associated Jacobian. The method presented 
in this paper can be used with any kind of reaction coor- 
dinate, cartesian coordinates or angles. Secondly, we do 
not have to guess the shape of the potential of mean force 
to guide the simulation with a biasing potential as in um- 
brella sampling. Instead, the method has its only input 
in the range of expected values of the reaction coordinate. 
It can subsequently provide the guiding potential which 
yields a uniform probability distribution. This "flat his- 
togram" allows one to get accurate estimates of various 
properties using reweighting techniques, especially when 
large free energy barriers are present. Thirdly, and in 
contrast with the adiabatic free energy dynamics of Rosso 
and coworkers |10, [HJ, no unknown parameter must be 
estimated for the method to work optimally. Fourthly, 
because it is for use with Monte Carlo simulations, it can 
take benefit of many ideas aiming at accelerating conver- 
gence, such as the preferential sampling scheme employed 
here for a solvated molecule. Lastly, but importantly, the 
method is not limited to any specific statistical ensem- 
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ble, since it requires only the knowledge of the equilib- 
rium statistical distribution. The very limited input of 
this algorithm makes it suitable for complex atomic and 
molecular systems, for which defining a suitable reaction 
coordinate can already be a difficult problem. 

We investigate three classes of model problems, namely 
a simple ten-dimensional double-well potential energy 
surface, the mean force of the dihedral angle in the 
gas phase and solvated n-butane molecule, and the 
cube^ring isomerisation of the rocksalt (NaF) 4 clus- 
ter. We obtain good agreement with previous works, 



when published data are available. The method can 
be straightforwardly extended to treat free-energy hyper- 
surfaces by considering multidimensional weighting func- 
tions g(Xi,...,X p ) and the corresponding histograms. 
This could also help in reaching equilibrium along co- 
ordinates, which are normal to the reaction coordinate. 
It is also a simple matter of algebra to generalise it to 
quantum free energy profiles by replacing the conven- 
tional Boltzmann factors by Feynman path integrals in 
the framework of quantum Monte Carlo. 
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